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Abstract.  We  propose  and  analyze  an  extremely  fast,  efficient  and  simple  method  for  solving  the  problem: 

min{||tx||i :  Au  =  f,u£  Rn}. 

This  method  was  first  described  in  [1],  with  more  details  in  [2]  and  rigorous  theory  given  in  [3]  and  [4].  The  motivation  was 
compressive  sensing,  which  now  has  a  vast  and  exciting  history,  which  seems  to  have  started  with  Candes,  et.al.  [5]  and  Donoho, 
[6].  See  [2],  [3]  and  [4]  for  a  large  set  of  references.  Our  method  introduces  an  improvement  called  “kicking”  of  the  very  efficient 
method  of  [1],  [2]  and  also  applies  it  to  the  problem  of  denoising  of  undersampled  signals.  The  use  of  Bregman  iteration  for 
denoising  of  images  began  in  [7]  and  led  to  improved  results  for  total  variation  based  methods.  Here  we  apply  it  to  denoise 
signals,  especially  essentially  sparse  signals,  which  might  even  be  undersampled. 
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1.  Introduction 

Let  A&Rmxn,  with  n>m  and  f£Rm,  be  given.  The  aim  of  a  basis  pursuit  problem  is  to  find  u&Rn 
by  solving  the  constrained  minimization  problem: 

min  {  J(u)\Au=  /}  (1.1) 

uS-R" 

where  J{u)  is  a  continuous  convex  function. 

For  basis  pursuit,  we  take: 

n 

J(u)  =  \u\i  =  ^2\uj\.  (1.2) 

i= i 

We  assume  that  AAT  is  invertible.  Thus  Au  =  f  is  underdetermined  and  has  at  least  one  solution,  u  = 
AT (AAT)~1  f ,  which  minimizes  the  ii  norm.  We  also  assume  that  J(u )  is  coercive,  i.e.,  whenever  ||m||— > 
oo,  J{u )  — >oo.  This  implies  that  the  set  of  all  solutions  of  (1.1)  is  nonempty  and  convex.  Finally,  when  J{u) 
is  strictly  or  strongly  convex,  the  solution  of  (1.1)  is  unique. 

Basis  pursuit  arises  from  many  applications.  In  particular,  there  has  been  a  recent  burst  of  research 
in  compressive  sensing,  which  involves  solving  (1.1),  (1.2).  This  was  led  by  Candes  et.al.  [5],  Donoho,  [6], 
and  others,  see  [2],  [3]  and  [4]  for  extensive  references.  Compressive  sensing  guarantees,  under  appropriate 
circumstances,  that  the  solution  to  (1.1),  (1-2)  gives  the  sparsest  solution  satisfying  Au=f.  The  problem 
then  becomes  one  of  solving  (1.1),  (1-2)  fast.  Conventional  linear  programming  solvers  are  not  tailored  for 
the  large  scale  dense  matrices  A  and  the  sparse  solutions  u  that  arise  here.  To  overcome  this,  a  linearized 
Bregman  iterative  procedure  was  proposed  in  [1]  and  analyzed  in  [2],  [3]  and  [4].  In  [2],  true,  nonlinear 
Bregman  iteration  was  also  used  quite  successfully  for  this  problem. 

Bregman  iteration  applied  to  (1.1),  (1-2)  involves  solving  the  constrained  optimization  problem  through 
solving  a  small  number  of  unconstrained  optimization  problems: 

min|/r|w|i  +  ^||Au-/||||  (1.3) 
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for  /i  >  0. 

In  [2]  we  used  a  method  called  the  fast  fixed  point  continuation  solver  (FPC)  [8]  which  appears  to  be 
efficient.  Other  solvers  of  (1.3)  could  be  used  in  this  Bregman  iterative  regularization  procedure. 

Here  we  will  improve  and  analyze  a  linearized  Bregman  iterative  regularization  procedure,  which,  in  its 
original  incarnation,  [1],  [2],  involved  only  a  two  line  code  and  simple  operations  and  was  already  extremely 
fast  and  accurate. 

In  addition,  we  are  interested  in  the  denoising  properties  of  Bregman  iterative  regularization,  for  signals, 
not  images.  The  results  for  images  involved  the  BV  norm,  which  we  may  discretize  for  n  x  n  pixel  images  as: 

n—  1 

TV (u)  =  ^2  ((ui+ltj-Uij)2  +  (uitj+1-Uij)2)K  (1.4) 

i,j= 1 

We  usually  regard  the  success  of  the  ROF  TV  based  model  [9] 

min{TV(U)  +  ^||/-u||2}  (1.5) 

(we  now  drop  the  subscript  2  for  the  L2  norm  throughout  the  paper)  as  due  to  the  fact  that  images  have 
edges  and  in  fact  are  almost  piecewise  constant  (with  texture  added).  Therefore,  it  is  not  surprising  that 
sparse  signals  could  be  denoised  using  (1.3).  The  ROF  denoising  model  was  greatly  improved  in  [7]  and  [10] 
with  the  help  of  Bregman  iterative  regularization.  We  will  do  the  same  thing  here  using  Bregman  iteration 
with  (1.3)  to  denoise  sparse  signals,  with  the  added  touch  of  undersampling  the  noisy  signals. 

The  paper  is  organized  as  follows:  In  section  2  we  describe  Bregman  iterative  algorithms,  as  well  as 
the  linearized  version.  We  motivate  these  methods  and  describe  previously  obtained  theoretical  results. 
In  section  3  we  introduce  an  improvement  to  the  linearized  version,  call  “kicking”  which  greatly  speeds 
up  the  method,  especially  for  solutions  u  with  a  large  dynamic  range.  In  section  4  we  present  numerical 
results,  including  sparse  recovery  for  u  having  large  dynamic  range,  and  the  recovery  of  signals  in  large 
amounts  of  noise.  In  another  work  in  progress  [11]  we  apply  these  ideas  to  denoising  very  blurry  and  noisy 
signals  remarkably  well  including  sparse  recovery  for  u.  By  blurry  we  mean  situations  where  A  is  perhaps  a 
subsampled  discrete  convolution  matrix  whose  elements  decay  to  zero  with  n,  e.g.  random  rows  of  a  discrete 
Gaussian. 

2.  Bregman  and  Linearized  Bregman  Iterative  Algorithms 

The  Bregman  distance  [12],  based  on  the  convex  function  J,  between  points  u  and  v,  is  defined  by 

Dj(u,v)  =  J(u)  —  J(v)  —  (p,u  —  v)  (2.1) 

where  p£dJ(v)  is  an  element  in  the  subgradient  of  J  at  the  point  v.  In  general  Dp(u,v)  ^  Dp(v,u)  and  the 
triangle  inequality  is  not  satisfied,  so  Dp(u,v)  is  not  a  distance  in  the  usual  sense.  However  it  does  measure 
the  closeness  between  u  and  v  in  the  sense  that  Dp(u,v )  >0  and  Dj(u,v )  >Dj(w,v )  for  all  points  w  on  the 
line  segment  connecting  u  and  v.  Moreover,  if  J  is  convex,  Dp(u,v)  >0,  if  J  is  strictly  convex  Dp(u,v)  >0 
for  u  7^  v  and  if  J  is  strongly  convex,  then  there  exists  a  constant  a  >  0  such  that 

Dj(u,v)  >a||u  — uj|2. 

To  solve  (1.1)  Bregman  iteration  was  proposed  in  [2]  .  Given  u°=p°  =  0,  we  define: 

ufc+1=argmin  {  J(u)  -  J(uk)  -  {u-uk ,pk) +  h\Au- f\\2}  (2.2) 

u£Rn  2  J 

pk+1  =pk  —  AT(Auk+1  —  /). 

This  can  be  written  as 

uk+1=  arg  min  {  Dp  (u,  uk)  +  \  \ \  Au  -  f  \ \ 2 1 . 
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It  was  proven  in  [2]  that  if  J(u)  £  C2(Cl)  and  is  strictly  convex  in  fi,  then  ||  Auk  —  /||  decays  exponentially 
whenever  uk  £  Cl  for  all  k.  Furthermore,  when  uk  converges,  its  limit  is  a  solution  of  (1-1)-  It  was  also  proven 
in  [2]  that  when  J(u)  =  i.e.  for  problem  (1.1)  and  (1.2),  or  when  J  is  a  convex  function  satisfying  some 
additional  conditions,  the  iteration  (2.2)  leads  to  a  solution  of  (1.1)  in  finitely  many  steps. 

As  shown  in  [2],  see  also  [7],  [10],  the  Bregman  iteration  (2.2)  can  be  written  as: 

fk+i  =  fk  +  f_Auk 

uk+1  =  arg  min  /  J(u)  +  ^-\\Au—  fk+1 1|21  (2.3) 

u€Rn  [  2  J 

This  was  referred  to  as  “adding  back  the  residual”  in  [7]  .  Here  f°  =  0,u°  =  0.  Thus  the  Bregman  iteration 
uses  solutions  of  the  unconstrained  problem 

min{^(w)  +  ^Pw-/H2}  (2-4) 

as  a  solver  in  which  the  Bregman  iteration  applies  this  process  iteratively. 

Since  there  is  generally  no  explicit  expression  for  the  solver  of  (2.2)  or  (2.3),  we  turn  to  iterative  methods. 
The  linearized  Bregman  iteration  which  we  will  analyze,  improve  and  use  here  is  generated  by 

uk+1  =  arg  min  (  J(u)  -  J(uk)  -  (u  -  uk ,pk)  +  J-  ||u ■ -  (uk  -  SAT ( Auk  -  /) )  || 2 1 
u&R71  (2  o  J 

pk+1  =  pk  _  i(ufc+1  -  uk )  -  AT(Auk  -  /).  (2.5) 

In  the  special  case  considered  here,  where  J(u)  =/z||u||i,  then  we  have  the  two  line  algorithm 

vk+1  =  vk~AT(Auk~f)  (2.6) 

uk+1  =  5  ■  shrink(ufc+1 ,  fi)  (2.7) 

where  vk  is  an  auxiliary  variable 

vk=pk+\  uk  (2.8) 

o 

and 

!x  —  fi,  if  x>  i-i 
0,  ii  —fi<x<n 
x  +  n,  if  x  <  —jjL 


is  the  soft  thresholding  algorithm  [13]  . 

This  linearized  Bregman  iterative  algorithm  was  invented  in  [1]  and  used  and  analyzed  in  [2], [3]  and  [4]. 
In  fact  it  comes  from  the  inner-outer  iteration  for  (2.2).  In  [2]  it  was  shown  that  the  linearized  Bregman 
iteration  (2.5)  is  just  one  step  of  the  inner  iteration  for  each  outer  iteration.  Here  we  repeat  the  arguments 
also  in  [2],  which  begin  by  summing  the  second  equation  in  (2.5)  arriving  at  (using  the  fact  that  u°  =p°  =  0): 

1  k~ 1  1 

pk  +  -uk  +  ^2AT(Auj  -  f)=pk  +  -uk -vk  =  0,  for  fc  =  l,2,.... 

3=0 

This  gives  us  (2.7),  and  allows  us  to  rewrite  its  first  equation  as: 

ufe!  1  =arg  min  /  J(u)+  v-j||w  —  5ufe+1||2l  (2.9) 

ueRn  [  2d  J 

i.e.  we  are  adding  back  the  “linearized  noise”,  where  vk+1  is  defined  in  (2.6). 
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In  [2]  and  [3]  some  interesting  analysis  was  done  for  (2.5),  (and  some  for  (2.9)).  This  was  done  first  for 
J(u)  continuously  differentiable  in  (2.5)  and  the  gradient  dJ(u)  satisfying 

||<9  J(u)  —  d  J(v)  1 1 2  <  P(d  J  (u)  —  dJ  (v)  ,u  —  v),  (2.10) 

\/u,vGRn\  /3>  0.  In  [3]  it  was  shown  that,  if  (2.10)  is  true,  then  both  of  the  sequences  ( uk)keN  and  ( pk)keN 
defined  by  (2.5)  converge  for  0  <  5  <  ■ 

In  [4]  the  authors  recently  give  a  theoretical  analysis,  showing  that  the  iteration  in  (2.6)  and  (2.7) 
converges  to  the  unique  solution  of 

min  {m||w||1  +  ^||w||2:-4m  =  /}  (2.11) 

They  also  show  the  interesting  result:  let  S  be  the  set  of  all  solutions  of  the  Basis  Pursuit  problem  (1.1), 
(1.2)  and  let 

u\  =  argmin||tt||2  (2.12) 

uGS 

which  is  unique.  Denote  the  solution  of  (2.11)  as  a*.  Then 

lim  ||u*  — Mill  =0.  (2.13) 

fl — >-00  ^ 

In  passing  they  show  that 

||w*||<||ui||  for  all  n>0  (2-14) 

which  we  will  use  below. 

Another  theoretical  analysis  on  Linearized  Bregman  algorithm  is  given  by  Yin  in  [14],  where  he  shows 
that  Linearized  Bregman  iteration  is  equivalent  to  gradient  descent  applied  to  the  dual  of  the  problem  (2.11) 
and  uses  this  fact  to  obtain  an  elegant  convergence  proof. 

This  summarizes  the  relevant  convergence  analysis  for  our  Bregman  and  linearized  Bregman  models. 
Next  we  recall  some  results  from  [7]  regarding  noise  and  Bregman  iteration. 

For  any  sequence  {uk},{pk}  satisfying  (2.2)  for  J  continuous  and  convex,  we  have,  for  any  p, 

D^(fi,Ufe)-DJpfc-1(w,Ufc-1)<(AfI-/,AUfc-1-/)-||AUM_/||2.  (2.15) 

Besides  implying  that  the  Bregman  distance  between  uk  and  any  element  u  satisfying  Au  =  /  is  mono- 
tonically  decreasing,  it  also  implies  that,  if  u  is  the  “noise  free”  approximation  to  the  solution  of  (1.1),  the 
Bregman  distance  between  uk  and  u  diminishes  as  long  as 

\\Auk  —  /||  >  ||  Aft—  /||  =  a,  where  a  is  some  measure  of  the  noise  (2.16) 

i.e.,  until  we  get  too  close  to  the  noisy  signal  in  the  sense  of  (2.16).  Note,  in  [7]  we  took  A  to  be  the  identity, 
but  these  more  general  results  are  also  proven  there.  This  gives  us  a  stopping  criterion  for  our  denoising 
algorithm. 

In  [7]  we  obtained  a  result  for  linearized  Bregman  iteration,  following  [15],  which  states  that  the  Bregman 
distance  between  u  and  uk  diminish  as  long  as 

\\Au-f\\<(l-2S\\AAT\\)  \\Auk  —  /||  (2.17) 

so  we  need  0  <  26||AAT||  <  1. 

In  practice,  we  will  use  (2.16)  as  our  stopping  criterion. 
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3.  Convergence 

We  begin  with  the  following  simple  results  for  the  linearized  Bregman  iteration  or  the  equivalent  algo¬ 
rithm  (2.5). 

Theorem  3.1.  Ifuk^u°°,  thenAu°°  =  f. 

Proof.  Assume  Au°°  yf  / .  Then  AT{Au°°  —  f)^  0  since  AT  has  full  rank.  This  means  that  for  some 
i,  (AT (Auk  —  /)),;  converges  to  a  nonzero  value,  which  means  that  vk+1  —  vk  does  as  well.  On  the  other 
hand  {vk}  =  {uk /5+pk}  is  bounded  since  {uk}  converges  and  pk  €  [— Therefore  {vk}  is  bounded,  while 
vk+1  —vk  converges  to  a  nonzero  limit,  which  is  impossible.  □ 

Theorem  3.2.  //wfc— >tt°°  andvk—>v°°,  then  u°°  minimizes  {J(u)+ ^\\u\\2  :  Au  =  f}. 

Proof.  Let  J(u)  =  J(u)  + ^\\u\\2.  then 

dJ(u)=dJ(u )  +  \u. 

o 

Since  dJ(uk)=pk=vk  —  uk/51  we  have  dJ(uk)  =  vk .  Using  the  non-negativity  of  the  Bregman  distance  we 
obtain 

J(uk)  <  j(uopt)  -  (uopt  -  uk,dJ{uk)) 

=  Jr('u0pt)  (^opt  ^  A  ) 

where  u0pt  minimizes  (1.1)  with  J  replaced  by  J,  which  is  strictly  convex. 

Let  A:  — >  oo,  we  have 

J (u°°)  <  J (uopt)  -  (Uopt  -  U°°  ,V°°) 

Since  vk  =  ATY^jZo  AT(/  —  Au^),  we  have  i>°°  Grange(AT).  Since  Auopt  =  Au°°  =  f ,  we  have  (uopt  — 
u°° ,v°°)  =  0,  which  implies  J(u°°)<J(uop t).  □  Equation  (2.11)  (from  a  result  in  [3]  )  implies  that  u°° 
will  approach  a  solution  to  (1.1),  (1.2),  as  fi  approaches  oo. 

The  linearized  Bregman  iteration  has  the  following  monotonicity  property: 

Theorem  3.3.  Ifuk+1^uk  and  0  <5  <2/\\AAT\\,  then 

\\Auk+1^f\\<\\Auk-f\\. 

Proof.  Let 

ufc+1-ufc  =  A«fc,  vk+1-vk  =  Avk. 

Then  the  shrinkage  operation  is  such  that 

Auk  =  5qkAvk  (3.1) 

with 

(  =  1  if  uk+luk  >  0 
9?  |=0  if  ui+1  =uk  =  0 
[e  (0,1]  otherwise 

Let  Qfc  =  Diag  (qk).  Then  (3.1)  can  be  written  as 

Auk  =  SQkAvk  =  SQkAT(f-Auk)  (3.2) 

which  implies 


Auk+1  -/=(/-  6AQk  At)  ( Auk  -  /) . 


(3.3) 
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From  (3.1),  Qk  is  diagonal  with  0  A  Qk  A  I,  so  0  ^  AQkAT  A  AAr .  If  we  choose  5  >  0  such  that  SAAT  A 
2 1,  then  0  A  5AQk  AT  A  2/  or  —I^I  —  5AQkATAI  which  implies  that  \\Auk  —  f\\  is  not  increasing.  To  get 
strict  decay,  we  need  only  show  that  AQk AT  (Auk  —  /)  =  0  is  impossible  if  uk+l  uk .  Suppose  AQkAT(Auk  — 
/)  =  0  holds,  then  from  (3.2)  we  have: 

(A  uk ,  Avk)  =  6(AT(f  -  Auk)  ,QkAT{f-Auk))  =  0. 

By  (3.1),  this  only  happens  if  uk+1  =uk  for  all  i,  which  is  a  contradiction.  □ 

We  are  still  faced  with  estimating  how  fast  the  residual  decays.  It  turns  out  that  if  consecutive  elements 
of  u  do  not  change  sign,  then  ||Au  —  f\\  decays  exponentially.  By  ’exponential’  we  mean  that  the  ratio  of  the 
residuals  of  two  consecutive  iteration  converges  to  a  constant,  this  type  of  convergence  is  sometimes  called 
linear  convergence.  Here  we  define 


S«  =  {a:e  Rn  ■  signal)  =  sign(uj) ,Vi}  (3.4) 

(where  sign(0)=0  and  sign(a)  =  a/\a\  for  ay^O).  Then  we  have  the  following: 

Theorem  3.4.  IfukGS  =  SUk  for  k£  {T\,T2),  thenvf  converges  ton*,  where  u*  €  argmin{||Hu  — /||2  :uGS} 
and  ||Aufe  — /||2  decays  to  ||Au*  — /||2  exponentially. 

Proof.  .  Since  uk  €  S  for  k  €  [Ti,T2],  we  can  define  Q  =  Qk  for  T\  <  k  <  T2  —  1.  From  (3.1)  we  see  that  Qk 
is  a  diagonal  matrix  consisting  of  zeros  or  ones,  so  Q  =  QTQ.  Moreover,  it  is  easy  to  see  that  S  =  {a’|Qx  =  x}. 
Following  the  argument  in  Theorem  3.3  we  have: 

uk+ 1  -uk  =  A  uk  =  6QAvk  =  6QAT(f  -  Auk )  (3.5) 

Auk+1-f=[I-SAQAT](Auk-f)  (3.6) 


and 


-I^I-SAQAt  <1- 

Let  Rn  =  V o  ©  Vi,  where  Vo  is  the  null  space  of  AQAT  and  V\  is  spanned  by  the  eigenvectors  corresponding 
to  the  nonzero  eigenvalues  of  AQAT.  Let  Auk  —  f  =  wk,° +  wk'1 ,  where  wk^  GVj  for  j  =  0,1.  From  (3.6)  we 
have 


wk+1’1  =  [I-5AQAT]wk’1 

for  Ti<k<T2  —  l.  Since  wk'1  is  not  in  the  null  space  of  AQAT,  then  (3.5)  and  (3.6)  imply  that  ||wfe,1| 
decays  exponentially.  Let  w°  =  wk,°,  then  AQAT w°  =  0  AQQATw°  =>QATw°  =  0.  Therefore,  from  (3.5)  we 
have 

A uk  =  6QTAT{f  -Auk)  =  6QAT(w°  +  wk -1)  =  SQATwk<1. 

Thus  ||  Aufe||  decays  exponentially.  This  means  {uk}  forms  a  Cauchy  sequence  in  S,  so  it  has  a  limit  u*  £  S. 
Moreover 


Au*  —  /  =  lim(Aufe  -  /)  =  lim  wk’°  +  lim  w^1  =  w°. 


Since  Vo  and  Vi  are  orthogonal: 

— /l|2  =  ||wfc,0||2  +  ||w/c,1||2  =  ||^4u*  —  /||2  +  ||wfc,1||2, 

so  \\Auk  —  f\\2  ~\\Au*  —  f\\2  decays  exponentially.  The  only  thing  left  to  show  is  that 

u*  =  argmin(||Hu  —  /||2  :u£  S)  =  argmin{||Hu  —  /  ||2  :  Qu=  u}. 

This  is  equivalent  to  way  that  AT(Au*  —  f)  is  orthogonal  with  the  hyperspace  {u:Qu  =  u}.  It’s  easy  to  see 
that  since  Q  is  a  projection  operator,  a  vector  v  is  orthogonal  with  {u :  Qu  =  it}  if  and  only  if  Qv  =  0,  thus  we 
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need  to  show  QAT(Au*  —  f)  =  0.  This  is  obvious  because  we  have  shown  that  Au*  —  f  =  w°  and  QATw°  =  0. 
So  we  find  that  u*  is  the  desired  minimizer.  □ 

Therefore,  instead  of  decaying  exponentially  with  a  global  rate,  the  residual  of  the  linearized  Bregman 
iteration  decays  in  a  rather  sophisticated  manner.  From  the  definition  of  the  shrinkage  function  we  can  see 
that  the  sign  of  an  element  of  u  will  change  if  and  only  if  the  corresponding  element  of  v  crosses  the  boundary 
of  the  interval  [—//,//].  If  /i  is  relatively  large  compared  with  the  size  of  Au  (which  is  usually  the  case  when 
applying  the  algorithm  to  a  compressed  sensing  problem),  then  at  most  iterations  the  signs  of  the  elements 
of  u  will  stay  unchanged,  i.e.  u  will  stay  in  the  subspace  Su  defined  in  (3.4)  for  a  long  while.  This  theorem 
tells  us  that  under  this  scenario  u  will  quickly  converge  to  the  point  u*  that  minimizes  ||Ait  — /||  inside  Su, 
and  the  difference  between  ||Au  —  f\\  and  ||Au*  —  /||  decays  exponentially.  After  u  converges  to  u* ,  it  will 
stay  there  until  the  sign  of  some  element  of  it  changes.  Usually  this  means  that  a  new  nonzero  element  of  it 
comes  up.  After  that,  it  will  enter  a  different  subspace  S  and  a  new  converging  procedure  begins. 

The  phenomenon  described  above  can  be  observed  clearly  in  Fig  3.1.  The  final  solution  of  u  contains 
five  non-zero  spikes.  Each  time  a  new  spike  appears,  it  converges  rapidly  to  the  position  that  minimizes 
||  Ait  — /||  in  the  subspace  Su.  After  that  there  is  a  long  stagnation,  which  means  it  is  just  waiting  there  until 
the  accumulating  v  brings  out  a  new  non-zero  element  of  u.  The  larger  /x  is,  the  longer  the  stagnation  takes. 
Although  the  convergence  of  the  residual  during  each  phase  is  fast,  the  total  speed  of  the  convergence  suffers 
much  from  the  stagnation.  The  solution  of  this  problem  will  be  described  in  the  next  section. 


0.8  - 

0.6  - 


10  20 


Fig.  3.1.  The  left  figure  presents  a  simple  signal  with  5  non-zero  spikes.  The  right  figure  shows  how  the  linearized  Bregman 
iteration  converges. 


4.  Fast  Implementation 

The  iterative  formula  in  Algorithm  1  below  gives  us  the  basic  linearized  Bregman  algorithm  designed  to 
solve  (1.1), (1.2). 


Algorithm  1  Bregman  Iterative  Regularization 
Initialize:  it  =  0,  v  =  0. 
while  “(I/  — Ait ||  not  converge”  do 
vk+1  =  vk  +  AT  (/  —  Auk) 
uk+1  =  <Tshrink(itfe+1,/x) 
end  while 


This  is  an  extremely  concise  algorithm,  simple  to  program,  involve  only  matrix  multiplication  and 
shrinkage.  When  A  consists  of  rows  of  a  matrix  of  a  fast  transform  like  FFT  which  is  a  common  case  for 
compressed  sensing,  it  is  even  faster  because  matrix  multiplication  can  be  implemented  efficiently  using  the 
existing  fast  code  of  the  transform.  Also,  storage  becomes  a  less  serious  issue. 


Fast  Linearized  Bregman  Iteration  for  Compressive  Sensing  and  Sparse  Denoising 


We  now  consider  how  we  can  accelerate  the  algorithm  under  the  problem  of  stagnation  described  in 
the  previous  section.  From  that  discussion,  during  a  stagnation  u  converges  to  a  limit  u*  so  we  will  have 
uk+1  uk+2  «  •  •  •  w uk+m  k, u*  for  some  m.  Therefore  the  increment  of  v  in  each  step,  AT(f  —  Au),  is  fixed. 
This  implies  that  during  the  stagnation  u  and  v  can  be  calculated  explicitly  as  following 


Uk+J=uk+1 

vk+i  =  vk  +  j  ■  Ar  (/  —  Auk+1) 


(4.1) 


If  we  denote  the  set  of  indices  of  the  zero  elements  of  u*  as  J0  and  let  I\  =  J0  be  the  support  of  u* ,  then  vk 
will  keep  changing  only  for  ie  Iq  and  the  iteration  can  be  formulated  entry- wise  as: 


fc+j  . 


Vi 


V-  =v!?+j-(Ar(f-Auk+1))i  iel0 

i  G  I\ 


k+j _  fc+1 

Va  =  Va 


(4.2) 


for  j  =  I,--  -  ,m.  The  stagnation  will  end  when  u  begins  to  change  again.  This  happens  if  and  only  if  some 
element  of  v  in  Iq  (which  keeps  changing  during  the  stagnation)  crosses  the  boundary  of  the  interval 
When  zG  J0,  vk  €  [— /z,/z],  so  we  can  estimate  the  number  of  the  steps  needed  for  vk  to  cross  the  boundary 
Vie  Jo  from  (4.2),  which  is 


sign ((AT(f  -  Auk+1))i )  -  vk+l 
(AT  (f  —  Auk+1))i 


Vie  J0 


(4.3) 


and 


s  =  min{sj} 
ie/o 


(4.4) 


is  the  number  of  steps  needed.  Therefore,  s  is  nothing  but  the  length  of  the  stagnation.  Using  (4.1),  we  can 
predict  the  end  status  of  the  stagnation  by 


uk+s  =  uk+ 1 

vk+s  =  vk+s-AT(f-Auk+1) 


(4.5) 


Therefore,  we  can  kick  u  to  the  critical  point  of  the  stagnation  when  we  detect  that  u  has  been  staying 
unchanged  for  a  while.  Specifically,  we  have  the  following  algorithm:  Algorithm  2. 


Algorithm  2  Linearized  Bregman  Iteration  with  Kicking 
Initialize:  u  =  0,  v  =  0. 
while  “||/  —  Au\\  not  converge”  do 
if  “ufc-1«  ufe”  then 

calculate  s  from  (4.3)  and  (4.4) 
vk+1  =  vk  +  s  ■  (AT (/  -  Auk))i,  Vi  e  J0 
vi=vh  Vie  Ji 
else 

vk+1  =  vk  +  AT(f-Auk) 

end  if 

uk+ 1  =  c)-shrink(t>fc+1,/i) 

end  while 


Indeed,  this  kicking  procedure  is  similar  to  line  search  commonly  used  in  optimization  problems  and 
modifies  the  initial  algorithm  in  no  way  but  just  accelerates  the  speed.  More  precisely,  note  that  the 
output  sequence  {uk,vk}  is  a  subsequence  of  the  original  one,  so  all  the  previous  theoretical  conclusions  on 
convergence  still  hold  here. 

An  example  of  the  algorithm  is  shown  in  Fig  4.1.  It  is  clear  that  all  the  stagnation  in  the  original 
convergence  collapses  to  single  steps.  The  total  amount  of  computation  is  reduced  dramatically. 
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Fig.  4.1.  The  left  figure  presents  the  convergence  curve  of  the  original  linearized  Bregman  iteration  using  the  same  signal 
as  Fig  3. 1 .  The  right  figure  shows  the  convergence  curve  of  the  linearized  Bregman  iteration  with  the  kicking  modification. 


5.  Numerical  Results  In  this  section,  we  demonstrate  the  effectiveness  of  the  algorithm  (with 
kicking)  in  solving  basis  pursuit  and  some  related  problems. 

5.1.  Efficiency 

Consider  the  constrained  minimization  problem 


min|u|i  s.t.  Au  =  f , 


where  the  constraints  Au  =  f  are  under-determined  linear  equations  with  A  an  to  x  n  matrix,  and  /  generated 
from  a  sparse  signal  u  that  has  a  number  of  nonzeros  re  <  to. 

Our  numerical  experiments  use  two  types  of  A  matrices:  Gaussian  matrices  whose  elements  were  gen¬ 
erated  from  i.i.d.  normal  distributions  7V"(0,1)  (randn(m,n)  in  MATLAB),  and  partial  discrete  cosine 
transform  (DCT)  matrices  whose  k  rows  were  chosen  randomly  from  the  nxn  DCT  matrix.  These  matrices 
are  known  to  be  efficient  for  compressed  sensing.  The  number  of  rows  to  is  chosen  as  to  ~  relog  (n/re)  for 
Gaussian  matrices  and  m~relogn  for  DCT  matrices  (following  [5]  ). 

The  tested  original  sparse  signals  u  had  numbers  of  nonzeros  equal  to  0.05 n  and  0.02n  rounded  to  the 
nearest  integers  in  two  sets  of  experiments,  which  were  obtained  by  round(0.05*n)  and  round(0.02*n)  in 
MATLAB,  respectively.  Given  a  sparsity  ||m||o,  he.,  the  number  of  nonzeros,  an  original  sparse  signal  Sgt" 
was  generated  by  randomly  selecting  the  locations  of  ||u||o  nonzeros,  and  sampling  each  of  these  nonzero 
elements  from  U{— 1,1)  (2*(rand-0.5)  in  MATLAB).  Then,  /  was  computed  as  Au.  When  ||u||o  is  small 
enough,  we  expect  the  basis  pursuit  problem,  which  we  solved  using  our  fast  algorithm,  to  yield  a  solution 
u*  =u  from  the  inputs  A  and  /. 

Note  that  partial  DCT  matrices  are  implicitly  stored  fast  transforms  for  which  matrix-vector  multi¬ 
plications  in  the  forms  of  Ax  and  ATx  were  computed  by  the  MATLAB  commands  dct(x)  and  idct(x), 
respectively.  Therefore,  we  were  able  to  test  on  partial  DCT  matrices  of  much  larger  sizes  than  Gaussian 
matrices.  The  sizes  m-by-n  of  these  matrices  are  given  in  the  first  two  columns  of  Table  5.1. 

Our  code  was  written  in  MATLAB  and  was  run  on  a  Windows  PC  with  a  Intel(R)  Core(TM)  2  Duo 
2.0GHz  CPU  and  2GB  memory.  The  MATLAB  version  is  7.4. 

The  set  of  computational  results  given  in  Table  5.1  was  obtained  by  using  the  stopping  criterion 


Uuk-f  || 
ll/ll 


<10"5, 


(5.1) 


which  was  sufficient  to  give  a  small  error  ||wfc  —  u||/||u||.  Throughout  our  experiments  in  Table  5.1,  we  used 
p=l  to  ensure  the  correctness  of  the  results. 
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Table  5.1.  Experiment  results  using  10  random  instances  for  each  configuration  of  (m,n,  ||u||o),  with  nonzero  elements 
of  u  come  from  U(— 1,1) . 


Results  of  linearized  Bregman-Li  with  kicking 


Stopping  tolerance. 

||  Auk 

-/ll/ll/ll<io-b 

Gaussian  matrices 

n 

m 

stopping  itr.  k 

relative  error  H^- 

-«ll/ll«ll 

time  (sec.) 

mean 

std. 

max 

mean 

std. 

max 

mean 

std. 

max 

||tt||o  =  0.05n 

1000 

300 

422 

67 

546 

2.0e-05 

4.3e-06 

2.7e-05 

0.42 

0.06 

0.51 

2000 

600 

525 

57 

612 

1.8e-05 

1.9e-06 

2.1e-05 

4.02 

0.45 

4.72 

4000 

1200 

847 

91 

1058 

1.7e-05 

1.7e-06 

1.9e-05 

25.7 

2.87 

32.1 

n 

m 

||u||o  =  0.02n 

1000 

156 

452 

98 

607 

2.3e-05 

2.6e-06 

2.6e-05 

0.24 

0.06 

0.33 

2000 

312 

377 

91 

602 

2.0e-05 

4.0e-06 

2.9e-05 

1.45 

0.38 

2.37 

4000 

468 

426 

30 

477 

1.6e-05 

2.1e-06 

2.0e-05 

6.96 

0.51 

7.94 

Partial  DCT  matrices 

n 

m 

||w||o  =  0.05n 

4000 

2000 

71 

6.6 

82 

9.1e-06 

2.5e-06 

1.2e-05 

0.43 

0.06 

0.56 

20000 

10000 

158 

14.5 

186 

6.2e-06 

2.1e-06 

l.le-05 

3.95 

0.36 

4.73 

50000 

25000 

276 

14 

296 

6.8e-06 

2.6e-06 

1.0e-05 

17.6 

0.99 

19.2 

n 

m 

|«||o  =  0.02n 

4000 

1327 

52 

7.0 

64 

8.6e-06 

1.3e-06 

l.le-05 

0.27 

0.04 

0.35 

20000 

7923 

91 

10.3 

115 

7.2e-06 

2.2e-06 

l.le-05 

2.36 

0.30 

3.02 

50000 

21640 

140 

9.7 

153 

5.9e-06 

2.4e-06 

l.le-05 

8.53 

0.66 

9.42 

5.2.  Robustness  to  Noise 

In  real  applications,  the  measurement  /  we  obtain  is  usually  contaminated  by  noise.  The  measurement 
we  have  is: 

/  =  /  +  n  =  Au  +  n,  nS7V(0,cr). 

To  characterize  the  noise  level,  we  shall  use  SNR.  (signal  to  noise  ratio)  instead  of  a  itself.  The  SNR  is 
defined  as  follows 

SNR(u)  :=201og10(jj^jj). 

In  this  section  we  test  our  algorithm  on  recovering  the  true  signal  u  from  A  and  the  noisy  measurement 
/.  As  in  the  last  section,  the  nonzero  entries  of  u  are  generated  from  IA{— 1,1),  and  A  is  either  a  Gaussian 
random  matrix  or  a  partial  DCT  matrix.  Our  stopping  criteria  is  given  by 

std (Auk  —  f)<a,  and  Iter.  <  1000, 

i.e.  we  stop  whenever  the  standard  deviation  of  residual  Auk  —  f  is  less  than  a  or  the  number  of  iterations 
exceeds  1000.  Table  5.2  shows  numerical  results  for  different  noise  level,  size  of  A  and  sparsity.  We  also  show 
one  typical  result  for  a  partial  DCT  matrix  with  size  n  =  4000  and  ||u||o  =  0.02n  =  80  in  Figure  5.1. 

5.3.  Recovery  of  Signal  with  High  Dynamical  Range  In  this  section,  we  test  our  algorithm 

on  signals  with  high  dynamical  ranges.  Precisely  speaking,  let  MAX  =  max{|wj|  :i  =  l,...,n}  and  MIN  = 
min{|wj|  =  The  signals  we  shall  consider  here  satisfy  ^jj^p«1010.  C*ur  u  is  generated  by 

multiplying  a  random  number  in  [0,1]  with  another  one  randomly  picked  from  {1,10,...,  1010}.  Here  we 
adopt  the  stopping  criteria 


<KTn 


\\Auk  —  /|| 
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Table  5.2.  Experiment  results  using  10  random  instances  for  each  configuration  of  (in,n,  ||u||o)  • 


Results  of  linearized  Bregman-Li  with  kicking 


Stopping  criteria. 

std(Awfc  -  /)  <  cr. 

Gaussian  matrices 

Avg.  SNR 

(n,m) 

stopping  itr.  k 

relative 

error  ||it* 

-sll/H 

time  (sec 

) 

mean 

stcl. 

max 

mean 

std. 

max 

mean 

std. 

max 

u||o  =  0.05n 

26.12 

(1000,300) 

420 

95 

604 

0.0608 

0.0138 

0.0912 

0.33 

0.09 

0.53 

25.44 

(2000,600) 

206 

32 

253 

0.0636 

0.0128 

0.0896 

1.49 

0.22 

1.79 

26.02 

(4000,1200) 

114 

11 

132 

0.0622 

0.0079 

0.0738 

3.32 

0.31 

3.81 

Avg.  SNR 

(n,m) 

u||o  =  0.02n 

27.48 

(1000,156) 

890 

369 

1612 

0.0456 

0.0085 

0.0599 

0.42 

0.17 

0.73 

25.06 

(2000,312) 

404 

64 

510 

0.0638 

0.0133 

0.0843 

1.37 

0.23 

1.74 

26.04 

(4000,468) 

216 

35 

267 

0.0557 

0.0068 

0.0639 

3.29 

0.55 

4.13 

Partial  DCT  matrices 

Avg.  SNR 

(n,m) 

it||o  =  0.05n 

23.97 

(4000,  2000) 

151 

9.2 

170 

0.0300 

0.0028 

0.0332 

0.94 

0.07 

1.03 

24.00 

(20000,10000) 

250 

14 

270 

0.0300 

0.0010 

0.0318 

7.88 

0.62 

8.86 

24.09 

(50000,25000) 

274 

9.9 

295 

0.0304 

0.0082 

0.0315 

20.4 

0.74 

20.1 

Avg.  SNR 

(n,m) 

u||o  =  0.02n 

24.29 

(4000,1327) 

130 

11 

157 

0.0223 

0.0023 

0.0253 

0.79 

0.08 

1.00 

24.37 

(20000,7923) 

223 

14 

257 

0.0204 

0.0025 

0.0242 

6.89 

0.53 

8.15 

24.16 

(50000,21640) 

283 

19 

311 

0.0193 

0.0012 

0.0207 

21.5 

1.68 

24.1 

Noisy  Measurements,  SNR  =  23.1084 

0. 

0. 

0. 

-0. 

-0. 

0  500  1000 


Iterations  =  102,  Final  Error  =  0.020764 


Fig.  5.1.  The  left  figure  presents  the  clean  (red  dots)  and  noisy  (blue  circles)  measurements,  with  SNR=23.108);  the 
right  figure  shows  the  reconstructed  signal  (blue  circles)  v.s.  original  signal  (red  dots),  where  the  relative  error=0. 02076),  and 
number  of  iterations  is  1 02. 


for  the  case  without  noise  (Figure  5.2)  and  the  same  stopping  criteria  as  in  the  previous  section  for  the  noisy 
cases  (Figures  5. 3-5. 5).  In  the  experiments,  we  take  the  dimension  n  =  4000,  the  number  of  nonzeros  of  u  to 
be  0.02n,  and  /t=1010.  Here  /i  is  chosen  to  be  much  larger  than  before,  because  the  dynamical  range  of  u 
is  large.  Figure  5.2  shows  results  for  the  noise  free  case,  where  the  algorithm  converges  to  a  10-11  residual 
in  less  than  300  iterations.  Figures  5. 3-5. 5  show  the  cases  with  noise  (the  noise  is  added  the  same  way  as  in 
previous  section).  As  one  can  see,  if  the  measurements  are  contaminated  with  less  noise,  signals  with  smaller 
magnitudes  will  be  recovered  well.  For  example  in  Figure  5.3,  the  SNR«118,  and  the  entries  of  magnitudes 
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104  are  well  recovered;  in  Figure  5.4,  the  SNR«97,  and  the  entries  of  magnitudes  105  are  well  recovered; 
and  in  Figure  5.5,  the  SNRa;49,  and  the  entries  of  magnitudes  107  are  well  recovered. 


x  10 


s  lter=  298. 


Zoom-In 
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Decay  of  log10(l|u-ut  Jl/Nu,  J 


Fig.  5.2.  Upper  left ,  true  signal  (red  dots)  v.s.  recovered  signal  (blue  circle);  upper  right ,  one  zoom-in  to  the  lower 
magnitudes;  lower  left,  decay  of  residual  log10  4  ;  lower  right,  decay  of  error  to  true  solution  log10  '"'l,  ;  ,u '  . 


5.4.  Recovery  of  Sinusoidal  Waves  in  Huge  Noise 

In  this  section  we  consider 


u(t)  =asin(at)  +  bcos(/3t)1 

where  a,b,a  and  (3  are  unknown.  The  observed  signal  u  is  noisy  and  has  the  form  u  =  u  +  n  with  n~Af(0,a). 
In  practice,  the  noise  in  u  could  be  huge,  i.e.  possibly  have  a  negative  SNR,  and  we  may  only  be  able 
to  observe  partial  information  of  u,  i.e.  only  a  subset  of  values  of  it  is  known.  Notice  that  the  signal  is 
sparse  (only  four  spikes)  in  frequency  domain.  Therefore,  this  is  essentially  a  compressed  sensing  problem 
and  fi-minimization  should  work  well  here.  Now  the  problem  can  be  stated  as  reconstructing  the  original 
signal  u  from  random  samples  of  the  observed  signal  u  using  our  fast  l-\  -minimization  algorithm.  In  our 
experiments,  the  magnitudes  a  and  b  are  generated  from  U{— 1,1);  frequencies  a  and  (3  are  random  multiples 
of  — ,  i.e.  a  =  k\  —  and  a  =  ki~,  with  kt  taken  from  {0,1,..., n  —  1}  randomly  and  n  denotes  the  dimension. 
We  let  I  be  a  random  subset  of  {1,2 ,...,?r}  and  f  =  u(I ),  and  take  A  and  AT  to  be  the  partial  matrix  of 
inverse  Fourier  matrix  and  Fourier  matrix  respectively.  Now  we  perform  our  algorithm  adopting  the  same 
stopping  criteria  as  in  section  5.2,  and  obtain  a  reconstructed  signal  denoted  as  x.  Notice  that  reconstructed 
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SNR  =  117.8948,  Iter  =  181,  Error  =  3.431 2e-007 
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Fig.  5.3.  Noisy  case.  Left  figure,  true  signal  (red  dots)  v.s.  recovered  signal  (blue  circle );  right  figure,  one  zoom-in  to  the 
magnitudes  105 .  The  error  is  measured  by  A*  . 


SNR  =  96.9765,  Iter  =  1 48,  Error  =  2.6386e-006 
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Fig.  5.4.  Noisy  case.  Left  figure,  true  signal  (red  dots)  v.s.  recovered  signal  (blue  circle);  right  figure,  one  zoom-in  to  the 
magnitudes  10® .  The  error  is  measured  by  . 


signal  x  is  in  Fourier  the  domain,  not  in  the  physical  domain.  Thus  we  take  an  inverse  Fourier  transform 
to  get  the  reconstructed  signal  in  physical  domain,  denoted  as  u* .  Since  we  know  a  priori  that  our  solution 
should  have  four  spikes  in  Fourier  domain,  before  we  take  the  inverse  Fourier  transform,  we  pick  the  four 
spikes  with  largest  magnitudes  and  set  the  rest  of  the  entries  to  be  zero.  Some  numerical  results  are  given  in 
Figure  7. 1-7.4.  Our  experiments  show  that  the  larger  the  noise  level  is,  the  more  random  samples  we  need  for 
a  reliable  reconstruction,  where  reliable  means  that  with  high  probability  (>80%)  of  getting  the  frequency 
back  exactly.  As  for  the  magnitudes  a  and  b ,  our  algorithm  cannot  guarantee  to  recover  them  exactly  (as 
one  can  see  in  Figure  7. 1-7. 4).  However,  frequency  information  is  much  more  important  than  magnitudes 
in  the  sense  that  the  reconstructed  signal  is  less  sensitive  to  errors  in  magnitudes  than  errors  in  frequencies 
(see  bottom  figures  in  Figure  7. 1-7. 4).  On  the  other  hand,  once  we  recover  the  right  frequencies,  one  can  use 
hardware  to  estimate  magnitudes  accurately. 
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SNR  =  49.289,  Iter  =  90,  Error  =  0.00067413 
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Fig.  5.5.  Noisy  case.  Left  figure,  true  signal  (red  dots)  v.s.  recovered  signal  (blue  circle);  right  figure,  one  zoom-in  to  the 
magnitudes  108 .  The  error  is  measured  by  _  7^  . 


6.  Conclusion  We  have  proposed  the  linearized  Bregman  iterative  algorithms  as  a  competitive 
method  for  solving  the  compressed  sensing  problem.  Besides  the  simplicity  of  the  algorithm,  the  special 
structure  of  the  iteration  enables  the  kicking  scheme  to  accelerate  the  algorithm  even  when  p  is  extremely 
large.  As  a  result,  a  sparse  solution  can  always  be  approached  efficiently. 

It  also  turns  out  that  our  process  has  remarkable  denoising  properties  for  undersampled  sparse  signals. 
We  will  pursue  this  in  further  work. 

Our  results  suggest  there  is  a  big  category  of  problem  that  can  be  solved  by  linearized  Bregman  iterative 
algorithms.  We  hope  that  our  method  and  its  extensions  could  produce  even  more  applications  for  problems 
under  different  scenarios,  including  very  underdetermined  inverse  problems  in  partial  differential  equations. 
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Orignal  and  Noisy  Waves,  SNR  =  2.6185 


Frequency  Domain 


Fig.  7.1.  Reconstruction  using  20%  random  samples  of  u  with  SNR=  2.6 1 85.  The  upper  left  figure  shows  the  original 
(red)  and  noisy  (blue)  signals;  the  upper  right  shows  the  reconstruction  (blue  circle)  v.s.  original  signal  (red  dots)  in  Fourier 
domain  in  terms  of  their  magnitudes  (i.e.  \u*\  v.s.  \u\);  bottom  left  shows  the  reconstructed  (blue)  v.s.  original  (red)  signal  in 
physical  domain;  and  bottom  right  shows  one  close-up  of  the  figure  at  bottom  left. 


[10]  T-C.  Chang,  L.  He,  and  T.  Fang.  Mr  image  reconstruction  from  sparse  radial  samples  using  bregman 
iteration.  Proceedings  of  the  13th  Annual  Meeting  of  ISMRM,  2006. 

[11]  Y.  Li,  S.  Osher,  and  Y.-H.  Tsai.  Recovery  of  sparse  noisy  date  from  solutions  to  the  heat  equation,  in 
preparation. 

[12]  L.M.  Bregman.  The  relaxation  method  of  finding  the  common  point  of  convex  sets  and  its  application  to 
the  solution  of  problems  in  convex  programming.  USSR  Computational  Mathematics  and  Mathematical 
Physics ,  7(3):200-217,  1967. 

[13]  D.L.  Donoho.  De-noising  by  soft-thresholding.  IEEE  Trans.  Inform.  Theory. 

[14]  W.  Yin.  On  the  linearized  bregman  algorithm,  private  communication. 

[15]  M.  Bachmayr.  Iterative  total  variation  methods  for  nonlinear  inverse  problems.  Master’s  thesis ,  Jo¬ 
hannes  Kepler  Universitat,  Linz,  Austria ,  2007. 


16 


Fast  Linearized  Bregman  Iteration  for  Compressive  Sensing  and  Sparse  Denoising 


Orignal  and  Noisy  Waves,  SNR  =  -4.7836 


Frequency  Domain 


Physical  Domain  One  Zoom-In 


Fig.  7.2.  Reconstruction  using  40%  random  samples  of  u  with  SNR=  —4.7836.  The  upper  left  figure  shows  the  original 
(red)  and  noisy  (blue)  signals;  the  upper  right  shows  the  reconstruction  (blue  circle)  v.s.  original  signal  (red  dots)  in  Fourier 
domain  in  terms  of  their  magnitudes  (i.e.  |u*|  v.s.  \u\);  bottom  left  shows  the  reconstructed  (blue)  v.s.  original  (red)  signal  in 
physical  domain;  and  bottom  right  shows  one  close-up  of  the  figure  at  bottom  left. 
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Orignal  and  Noisy  Waves,  SNR  =  -6.7908 
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Fig.  7.3.  Reconstruction  using  60%  random  samples  of  u  with  SNR=  —6.7908.  The  upper  left  figure  shows  the  original 
(red)  and  noisy  (blue)  signals;  the  upper  right  shows  the  reconstruction  (blue  circle)  v.s.  original  signal  (red  dots)  in  Fourier 
domain  in  terms  of  their  magnitudes  (i.e.  |u*|  v.s.  \u\);  bottom  left  shows  the  reconstructed  (blue)  v.s.  original  (red)  signal  in 
physical  domain;  and  bottom  right  shows  one  close-up  of  the  figure  at  bottom  left. 
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Orignal  and  Noisy  Waves,  SNR  =  -1 1 .001 6 
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Fig.  7.4.  Reconstruction  using  80%  random  samples  of  u  with  SNR=  —11.0016.  The  upper  left  figure  shows  the  original 
(red)  and  noisy  (blue)  signals;  the  upper  right  shows  the  reconstruction  (blue  circle)  v.s.  original  signal  (red  dots)  in  Fourier 
domain  in  terms  of  their  magnitudes  (i.e.  |n*|  v.s.  \u\);  bottom  left  shows  the  reconstructed  (blue)  v.s.  original  (red)  signal  in 
physical  domain;  and  bottom  right  shows  one  close-up  of  the  figure  at  bottom  left. 


